
import numpy as np
import pandas as pd
import geopandas as gpd
import fiona
import glob
import os
import contextily as ctx
from scipy.spatial import cKDTree
from shapely.geometry import Point, shape, LineString, MultiLineString, GeometryCollection, MultiPoint, Polygon, MultiPolygon  # creating points
import json
from tqdm.auto import tqdm
pd.set_option('min_rows', 30)
import sys
from importlib import reload
from util import *
import matplotlib.pyplot as plt
# from pandarallel import pandarallel # parallel operations for speed
# pandarallel.initialize(nb_workers=8, progress_bar=True)
import swifter
tqdm.pandas()
plt.rcParams['figure.figsize'] = (16, 16)
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 10)


ls('restricted')

Total: 0.0MB


%%time
parcels = gpd.read_file('input/lds-nz-primary-parcels-CLIPPED-4326.gpkg')

/data/miniconda3/envs/house-upzone/lib/python3.8/site-packages/geopandas/geodataframe.py:577: RuntimeWarning: Sequential read of iterator was interrupted. Resetting iterator. This can negatively impact the performance.
  for feature in features_lst:

CPU times: user 1min 22s, sys: 4.92 s, total: 1min 27s
Wall time: 1min 27s


parcels_sample = parcels.sample(10000)


# keep a centroid and polygon version - change between them according to need
parcels_sample['geometry_polygon_2193'] = parcels_sample.geometry.to_crs(2193)
parcels_sample['geometry_centroid_2193'] = parcels_sample.geometry_polygon_2193.centroid
parcels_sample['geometry_centroid_4326'] = parcels_sample.geometry_centroid_2193.to_crs(4326)
parcels_sample['geometry_polygon_4326'] = parcels_sample.geometry_polygon_2193.to_crs(4326)


bounds2193 = {'x1': 1.5e6, 'x2':2e6, 'y1':5.8e6, 'y2':6.1e6}


parcels_sample2 = parcels_sample[parcels_sample.parcel_intent != 'Road'].sample(100)


parcels_sample['LINZ_parcel_ID'] = parcels_sample.id


parcels_sample['geometry'] = parcels_sample.geometry_centroid_4326
parcels_sample['LINZ_parcel_centroid_lon'] = parcels_sample.geometry.x
parcels_sample['LINZ_parcel_centroid_lat'] = parcels_sample.geometry.y


parcels_sample['geometry'] = parcels_sample.geometry_polygon_4326
# all are multipolygons
print(set([type(a) for a in parcels_sample.geometry]))

{<class 'shapely.geometry.multipolygon.MultiPolygon'>}


%%time
# one way
coordinates_all = []
for multipolygon in tqdm(parcels_sample.geometry):
    coordinates = []
    for polygon in multipolygon:
        coordinates.extend(polygon.exterior.coords)
    coordinates_all.append(coordinates)

CPU times: user 1.44 s, sys: 28.1 ms, total: 1.47 s
Wall time: 1.46 s


coordinates_all[0]

[(174.67816820000007, -36.66715249999993),
 (174.6782238500001, -36.66718856699999),
 (174.67730413300015, -36.66768149999995),
 (174.67725030000008, -36.66764448299995),
 (174.67816820000007, -36.66715249999993)]


%%time
parcels_sample['LINZ_parcel_vertices_lon'] = [[p[0] for p in points] for points in coordinates_all]
parcels_sample['LINZ_parcel_vertices_lat'] = [[p[1] for p in points] for points in coordinates_all]

CPU times: user 115 ms, sys: 7.97 ms, total: 123 ms
Wall time: 122 ms


parcels_sample['LINZ_parcel_vertices_lon']

332148    [174.67816820000007, 174.6782238500001, 174.67...
138662    [174.74038025000004, 174.74054571700003, 174.7...
276820    [174.92074408300007, 174.92075735000003, 174.9...
31004     [174.70060606700008, 174.7007802500001, 174.70...
494255    [174.744322333, 174.7446516330001, 174.7445430...
                                ...                        
527367    [175.07140310000014, 175.07167341700006, 175.0...
24891     [174.49126546700006, 174.49126696700012, 174.4...
25636     [174.9080082170001, 174.90820755000004, 174.90...
369084    [174.7215236500001, 174.72155736700006, 174.72...
334475    [174.59673236700007, 174.5966783330001, 174.59...
Name: LINZ_parcel_vertices_lon, Length: 10000, dtype: object


parcels_sample[["LINZ_parcel_vertices_lon", "LINZ_parcel_vertices_lat"]].iloc[0].LINZ_parcel_vertices_lat

'[-36.667152; -36.667189; -36.667681; -36.667644; -36.667152]'


%%time
# another way (from phase 3)
def extract_verts(geometry):
    lat = np.nan
    lng = np.nan
    if geometry:
        coordinates = []
        for polygon in geometry:
            coordinates.extend(polygon.exterior.coords)
        lng = f"[{'; '.join([str(round(point[0], 6)) for point in coordinates])}]"
        lat = f"[{'; '.join([str(round(point[1], 6)) for point in coordinates])}]"
    return pd.Series([lng, lat])
parcels_sample[["LINZ_parcel_vertices_lon", "LINZ_parcel_vertices_lat"]] = parcels_sample.geometry.apply(extract_verts)

CPU times: user 6.9 s, sys: 43.3 ms, total: 6.94 s
Wall time: 6.93 s


roads = parcels[parcels.parcel_intent == "Road"]
roads = roads.to_crs(4326)


%%time
roads_dissolved = roads.dissolve()


parcels_sample2['geometry'] = parcels_sample2.geometry_polygon_4326
parcels_sample2 = parcels_sample2.set_crs(4326)


# faster than below
def get_points_in_roads(geom):
    """return a list of points from the geometry that fall within roads_dissolved"""
    # assume multipolygon
    assert isinstance(geom, MultiPolygon)
    assert len(geom) == 1
    for poly in geom:
        coords = poly.exterior.coords
        pointsx = [x for x, _ in coords]
        pointsy = [y for _, y in coords]
        points_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(pointsx, pointsy)).set_crs(4326)
        road_points = gpd.sjoin(points_gdf, roads_dissolved, op = 'intersects').geometry
    return road_points.geometry.values
#         return (road_points.x.values, road_points.y.values)

    # 
# within_points = [get_points_in_roads(g) for g in parcels_sample2.geometry]
parcels_sample2.geometry.progress_apply(get_points_in_roads)

353059    [POINT (175.054731883 -37.07730548299998), POI...
148404       [POINT (174.6879781000001 -36.80336151699999)]
506845       [POINT (174.7554596330001 -36.77824928299997)]
314555           [POINT (174.726871283 -36.92904104999997)]
180722    [POINT (174.596338583 -36.87381121699996), POI...
Name: geometry, dtype: object


%%time
def extract_road_intersections(geom):
    if geom:
        xmin, ymin, xmax, ymax = geom.bounds
        neighbouring_roads = roads.cx[xmin:xmax, ymin:ymax]
        return neighbouring_roads.geometry.intersection(geom).unary_union
    return np.nan
road_intersections = parcels_sample2.geometry_polygon_4326.progress_apply(extract_road_intersections)

CPU times: user 19.4 s, sys: 15.8 ms, total: 19.4 s
Wall time: 19.4 s


%%time
def extract_road_intersections2(geom):
    if geom:
        neighbouring_roads = roads.cx[xmin:xmax, ymin:ymax]
        return neighbouring_roads.geometry.intersection(geom).unary_union
    return np.nan
road_intersections = parcels_sample2.geometry_polygon_4326.progress_apply(extract_road_intersections)


types = pd.Series(type(i) for i in road_intersections)
types.index = road_intersections.index
types.value_counts()

<class 'shapely.geometry.collection.GeometryCollection'>      31
<class 'shapely.geometry.polygon.Polygon'>                    30
<class 'shapely.geometry.linestring.LineString'>              18
<class 'shapely.geometry.point.Point'>                        13
<class 'shapely.geometry.multipoint.MultiPoint'>               4
<class 'shapely.geometry.multipolygon.MultiPolygon'>           3
<class 'shapely.geometry.multilinestring.MultiLineString'>     1
dtype: int64


%%time
LINZ_parcel_roadvertices_lon = []
LINZ_parcel_roadvertices_lat = []
for intersection in tqdm(road_intersections):
    lat = np.nan
    lon = np.nan
    if intersection:
        lat = set()
        lon = set()
        try:
            if type(intersection) in (LineString, Point):
                lon.update(intersection.xy[0])
                lat.update(intersection.xy[1])
            elif type(intersection) in (MultiLineString, GeometryCollection, MultiPoint):
                for p in intersection.geoms:
                    if type(p) == Polygon:
                        lon.update(p.exterior.xy[0])
                        lat.update(p.exterior.xy[1])
                    else:
                        lon.update(p.xy[0])
                        lat.update(p.xy[1])
            else:
                raise Exception(f"Don't know how to handle {type(intersection)}")
        except:
            print(type(intersection))
            print(intersection)
            raise
        lon = f"[{'; '.join([str(round(lon, 6)) for lon in list(lon)])}]"
        lat = f"[{'; '.join([str(round(lat, 6)) for lat in list(lat)])}]"
    LINZ_parcel_roadvertices_lon.append(lon)
    LINZ_parcel_roadvertices_lat.append(lat)
# df["LINZ_parcel_roadvertices_lon"] = LINZ_parcel_roadvertices_lon
# df["LINZ_parcel_roadvertices_lat"] = LINZ_parcel_roadvertices_lat

<class 'shapely.geometry.polygon.Polygon'>
POLYGON ((174.7763490000001 -36.94952951699997, 174.7763300099195 -36.94948762124847, 174.7763490000001 -36.94952951699998, 174.7763490000001 -36.94952951699999, 174.7763490000001 -36.94952951699997))

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<timed exec> in <module>

Exception: Don't know how to handle <class 'shapely.geometry.polygon.Polygon'>


!ls restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/*power*.shp

ls: cannot access 'restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/*power*.shp': No such file or directory


power = gpd.read_file('input/Transmission_Lines_exTRANSPOWER/Transmission_Lines.shp')
power = power.cx[bounds2193['x1']:bounds2193['x2'], bounds2193['y1']:bounds2193['y2']]
# only interested in overhead
power = power[power['type'] == 'TRANSLINE']
power.sample(3)


power


ax = power.plot(column='designvolt', legend=True)
plt
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


parcels_sample['geometry'] = parcels_sample['geometry_polygon_2193']
parcels_sample = parcels_sample.set_crs(2193)
# get dataframe associating parcel indices with overhead power lines
parcels_sample_under = gpd.sjoin(parcels_sample, power[['descriptio', 'geometry']]).drop(columns=['index_right'])


parcels_sample_under.index.value_counts()

379134    3
373169    3
477337    3
332672    2
387854    2
         ..
444457    1
219029    1
57635     1
323105    1
228991    1
Length: 95, dtype: int64


id = 477337
ax = parcels_sample_under.loc[id].plot()
bounds = parcels_sample_under.loc[id].iloc[[0]].bounds
plt.xlim((bounds.minx.values[0] - 100, bounds.maxx.values[0] + 100))
plt.ylim((bounds.miny.values[0] - 100, bounds.maxy.values[0] + 100))
power.cx[bounds.minx.values[0]:bounds.maxx.values[0], bounds.miny.values[0]:bounds.maxy.values[0]].plot(column='descriptio', ax = ax, legend=True)

<AxesSubplot:>


!ls restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/*View[Ss]haf*.shp

restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_AucklandMuseumViewshaftContoursOverlay.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_AucklandMuseumViewshaftOverlay.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_DilworthTerraceHousesViewshaftContoursOverlay.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_DilworthTerraceHousesViewshaftOverlay.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_LocallySignificantVolcanicViewshaftsContours.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_LocallySignificantVolcanicViewshafts.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_RegionallySignificantVolcanicViewShaftsAndHeightSensitiveAreasOverlay.shp
restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_RegionallySignificantVolcanicViewshaftsContours.shp


!ls restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/*ViewShaf*.shp

restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_RegionallySignificantVolcanicViewShaftsAndHeightSensitiveAreasOverlay.shp


viewshafts_local = gpd.read_file('restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_LocallySignificantVolcanicViewshafts.shp').to_crs(2193)
viewshafts_regional = gpd.read_file('restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_RegionallySignificantVolcanicViewShaftsAndHeightSensitiveAreasOverlay.shp').to_crs(2193)


ax = viewshafts_regional[viewshafts_regional.NAME.isin(viewshafts_local.NAME.unique())].plot(color='blue', alpha=0.6)
viewshafts_local.plot(ax=ax, color='red', alpha=0.6)
blue_patch = mpatches.Patch(color='blue', label='regionally significant viewshafts')
red_patch = mpatches.Patch(color='red', label='locally significant viewshafts')
plt.legend(handles=[red_patch, blue_patch])
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


ax = viewshafts_local.plot(column='NAME', alpha=0.6, legend=True)
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


ax = viewshafts_regional.plot(column='NAME', alpha=0.6, legend=True)
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


viewshafts_museum = gpd.read_file('restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_AucklandMuseumViewshaftOverlay.shp').to_crs(2193)
viewshafts_dilworth = gpd.read_file('restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_DilworthTerraceHousesViewshaftOverlay.shp').to_crs(2193)


ax = viewshafts_museum.plot(alpha=0.6, legend=True)
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


ax = viewshafts_dilworth.plot(alpha=0.6, legend=True)
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=2193)


parcels_sample['geometry'] = parcels_sample.geometry_polygon_2193
parcels_sample = parcels_sample.set_crs(2193)
gpd.sjoin(parcels_sample, viewshafts_regional, how='left')


def extract_viewshaft_info(


parcels_sample.geometry.iloc[0].overlaps()

<bound method BaseGeometry.overlaps of <shapely.geometry.multipolygon.MultiPolygon object at 0x7f8c3ef45fa0>>


%%time
aup_zones = gpd.read_file('restricted/2016_Unitary_Plan_operational_15Nov/UP_OIP_15Nov2016_SHP/MASTER_UP_BaseZone.shp')
aup_zones = aup_zones.to_crs(2193)
aup_zones.sample(3)

CPU times: user 31.2 s, sys: 2.21 s, total: 33.4 s
Wall time: 33.3 s


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample = gpd.sjoin(parcels_sample, aup_zones[['ZONE', 'ZONE_resol', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'ZONE': 'LINZmatch_AUP_code', 'ZONE_resol': 'LINZmatch_AUP_name'})

CPU times: user 59.9 s, sys: 0 ns, total: 59.9 s
Wall time: 1min


ax = parcels_sample[parcels_sample.LINZmatch_AUP_name.str.contains('Business')].sample(500).plot(alpha=0.5)
ctx.add_basemap(ax, crs=2193)


rural_codes = aup_zones[aup_zones.ZONE_resol.str.lower().str.contains('rural - ', na=False)]['ZONE'].unique()


# check that each rural zone code matches with a unique rural zone name 
assert all([len(aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()) == 1 for code in rural_codes])
# dictionary mapping code to names
rural_code2name = {code: aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()[0] for code in rural_codes}


aup_zones[aup_zones.ZONE_resol.isna()]


# 2 NAs in ZONE_resol are from a zone 58, which only has observations
aup_zones[aup_zones.ZONE == '58']


rural = aup_zones[aup_zones.ZONE.isin(rural_codes)]


%%time
# method 1, slightly slower: dissolve by zone
rural_by_zone = rural.dissolve(by='ZONE').reset_index(drop=False)
distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for index, rural_row in rural_by_zone.iterrows():
        distance_candidates.append(row.geometry.distance(rural_row.geometry))
        code_candidates.append(rural_row.ZONE)
    distances.append(distance_candidates)
    codes.append(code_candidates)
    
# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 1min 2s, sys: 363 ms, total: 1min 3s
Wall time: 1min 2s


%%time
# method 2, slightly faster: dissolve all after subsetting gpd dfs
rural_by_zone_dict = {code: rural[rural.ZONE == code].dissolve() for code in rural_codes}

distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for code, rural_gdf in rural_by_zone_dict.items():
        distance_candidates.append(row.geometry.distance(rural_gdf.geometry[0]))
        code_candidates.append(rural_gdf.ZONE[0])
    distances.append(distance_candidates)
    codes.append(code_candidates)

# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 56.2 s, sys: 289 ms, total: 56.5 s
Wall time: 56.2 s


parcels_sample['Hdist_rural'] = distances[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_rural_code'] = codes[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_rural_name'] = parcels_sample.apply(lambda x: rural_code2name[x.Hdist_rural_code], axis=1)


import matplotlib.patches as mpatches
import matplotlib.pyplot as plt

colours = ('blue', 'orange', 'green', 'red', 'purple', 'brown', 'pink', 'cyan')
name2colour = {name: colour for name, colour in zip(rural_code2name.values(), colours)}

column = 'Hdist_rural'
subsample = parcels_sample[parcels_sample[column] > 0].sample(5)
subsample['buffer'] = subsample.apply(lambda x: x.geometry.buffer(x[column]), axis=1)
subsample['geometry'] = subsample['buffer']
subsample = subsample[~subsample.is_empty]
# ax = rural.plot(column='ZONE_resol', legend=True)
# subsample.plot(column='Hdist_rural_name', alpha=0.4, ax=ax)
ax = rural.plot(color=[name2colour[z] for z in rural.ZONE_resol])
subsample.plot(color=[name2colour[z] for z in subsample.Hdist_rural_name], alpha=0.65, ax=ax)
plt.ylim((5.88e6, 5.96e6))

# add legend
legend_patches = [mpatches.Patch(color=col, label=name) for name, col in name2colour.items()]
plt.legend(handles=legend_patches)

ctx.add_basemap(ax=ax, crs=2193)


import matplotlib
matplotlib.colors

<module 'matplotlib.colors' from '/data/miniconda3/envs/house-upzone/lib/python3.8/site-packages/matplotlib/colors.py'>


business_codes = aup_zones[aup_zones.ZONE_resol.str.lower().str.contains('business - ', na=False)]['ZONE'].unique()


# check that each business zone code matches with a unique business zone name 
assert all([len(aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()) == 1 for code in business_codes])
# dictionary mapping code to names
business_code2name = {code: aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()[0] for code in business_codes}
business_code2name

{'44': 'Business - Neighbourhood Centre Zone',
 '12': 'Business - Mixed Use Zone',
 '17': 'Business - Light Industry Zone',
 '5': 'Business - Heavy Industry Zone',
 '49': 'Business - General Business Zone',
 '1': 'Business - Business Park Zone',
 '22': 'Business - Town Centre Zone',
 '10': 'Business - Metropolitan Centre Zone',
 '7': 'Business - Local Centre Zone',
 '35': 'Business - City Centre Zone'}


business = aup_zones[aup_zones.ZONE.isin(business_codes)]


%%time
# method 1, slightly slower: dissolve by zone
business_by_zone = business.dissolve(by='ZONE').reset_index(drop=False)
distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for index, business_row in business_by_zone.iterrows():
        distance_candidates.append(row.geometry.distance(business_row.geometry))
        code_candidates.append(business_row.ZONE)
    distances.append(distance_candidates)
    codes.append(code_candidates)
    
# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 1min 2s, sys: 363 ms, total: 1min 3s
Wall time: 1min 2s


%%time
# method 2, slightly faster: dissolve all after subsetting gpd dfs
business_by_zone_dict = {code: business[business.ZONE == code].dissolve() for code in business_codes}

distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for code, business_gdf in business_by_zone_dict.items():
        distance_candidates.append(row.geometry.distance(business_gdf.geometry[0]))
        code_candidates.append(business_gdf.ZONE[0])
    distances.append(distance_candidates)
    codes.append(code_candidates)

# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 44.9 s, sys: 230 ms, total: 45.1 s
Wall time: 44.8 s


parcels_sample['Hdist_bus'] = distances[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_bus_code'] = codes[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_bus_name'] = parcels_sample.apply(lambda x: business_code2name[x.Hdist_bus_code], axis=1)


import matplotlib.patches as mpatches
import matplotlib.pyplot as plt

colours = ('blue', 'orange', 'green', 'red', 'purple', 'brown', 'pink', 'cyan', 'black', 'white')
name2colour = {name: colour for name, colour in zip(business_code2name.values(), colours)}

# hard to see, subset to smaller area
bounds = {'x1': 1.76e6, 'x2': 1.7613e6, 'y1': 5.9123e6, 'y2': 5.914e6}

column = 'Hdist_bus'
subsample = parcels_sample.cx[bounds['x1']:bounds['x2'], bounds['y1']:bounds['y2']]
subsample['buffer'] = subsample.apply(lambda x: x.geometry.buffer(x[column]), axis=1)
subsample['geometry'] = subsample['buffer']

subsample = subsample[~subsample.is_empty]

business_plot = business.cx[bounds['x1']:bounds['x2'], bounds['y1']:bounds['y2']]
ax = business_plot.plot(color=[name2colour[z] for z in business_plot.ZONE_resol])
subsample.plot(color=[name2colour[z] for z in subsample.Hdist_bus_name], alpha=0.65, ax=ax)
# plt.ylim((5.88e6, 5.96e6))

# add legend
legend_patches = [mpatches.Patch(color=col, label=name) for name, col in name2colour.items()]
plt.legend(handles=legend_patches)

ctx.add_basemap(ax=ax, crs=2193)

/data/miniconda3/envs/house-upzone/lib/python3.8/site-packages/geopandas/geodataframe.py:1322: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super(GeoDataFrame, self).__setitem__(key, value)


resid_codes = aup_zones[aup_zones.ZONE_resol.str.lower().str.contains('resid', na=False)]['ZONE'].unique()


resid_codes

array(['60', '23', '18', '8', '19', '20'], dtype=object)


# check that each resid zone code matches with a unique resid zone name 
assert all([len(aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()) == 1 for code in resid_codes])
# dictionary mapping code to names
resid_code2name = {code: aup_zones[aup_zones.ZONE == code].ZONE_resol.unique()[0] for code in resid_codes}
resid_code2name

{'60': 'Residential - Mixed Housing Urban Zone',
 '23': 'Residential - Large Lot Zone',
 '18': 'Residential - Mixed Housing Suburban Zone',
 '8': 'Residential - Terrace Housing and Apartment Building Zone',
 '19': 'Residential - Single House Zone',
 '20': 'Residential - Rural and Coastal Settlement Zone'}


resid = aup_zones[aup_zones.ZONE.isin(resid_codes)]


%%time
# method 1, slightly slower: dissolve by zone
resid_by_zone = resid.dissolve(by='ZONE').reset_index(drop=False)
distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for index, resid_row in resid_by_zone.iterrows():
        distance_candidates.append(row.geometry.distance(resid_row.geometry))
        code_candidates.append(resid_row.ZONE)
    distances.append(distance_candidates)
    codes.append(code_candidates)
    
# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 1min 2s, sys: 363 ms, total: 1min 3s
Wall time: 1min 2s


%%time
# method 2, slightly faster: dissolve all after subsetting gpd dfs
resid_by_zone_dict = {code: resid[resid.ZONE == code].dissolve() for code in resid_codes}

distances = []
codes = []
for index, row in tqdm(parcels_sample.iterrows(), total=len(parcels_sample)):
    distance_candidates = []
    code_candidates = []
    for code, resid_gdf in resid_by_zone_dict.items():
        distance_candidates.append(row.geometry.distance(resid_gdf.geometry[0]))
        code_candidates.append(resid_gdf.ZONE[0])
    distances.append(distance_candidates)
    codes.append(code_candidates)

# all distances (to any zone)
distances = np.array(distances)
codes = np.array(codes)
# indices of minimum distances
min_idx = np.argmin(distances, axis=-1)

CPU times: user 2min 9s, sys: 509 ms, total: 2min 9s
Wall time: 2min 8s


parcels_sample['Hdist_resid'] = distances[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_resid_code'] = codes[(np.arange(len(distances)), min_idx)]
parcels_sample['Hdist_resid_name'] = parcels_sample.apply(lambda x: resid_code2name[x.Hdist_resid_code], axis=1)


import matplotlib.patches as mpatches
import matplotlib.pyplot as plt

colours = ('blue', 'orange', 'green', 'red', 'purple', 'brown', 'pink', 'cyan', 'black', 'white')
name2colour = {name: colour for name, colour in zip(resid_code2name.values(), colours)}

# hard to see, subset to smaller area
bounds = {'x1': 1.76e6, 'x2': 1.7613e6, 'y1': 5.9123e6, 'y2': 5.914e6}

column = 'Hdist_resid'
subsample = parcels_sample.cx[bounds['x1']:bounds['x2'], bounds['y1']:bounds['y2']]
subsample['buffer'] = subsample.apply(lambda x: x.geometry.buffer(x[column]), axis=1)
subsample['geometry'] = subsample['buffer']

subsample = subsample[~subsample.is_empty]

resid_plot = resid.cx[bounds['x1']:bounds['x2'], bounds['y1']:bounds['y2']]
ax = resid_plot.plot(color=[name2colour[z] for z in resid_plot.ZONE_resol])
subsample.plot(color=[name2colour[z] for z in subsample.Hdist_resid_name], alpha=0.65, ax=ax)
# plt.ylim((5.88e6, 5.96e6))

# add legend
legend_patches = [mpatches.Patch(color=col, label=name) for name, col in name2colour.items()]
plt.legend(handles=legend_patches)

ctx.add_basemap(ax=ax, crs=2193)

/data/miniconda3/envs/house-upzone/lib/python3.8/site-packages/geopandas/geodataframe.py:1322: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super(GeoDataFrame, self).__setitem__(key, value)


%%time
postfix2name = {
            'SH': 'Residential - Single House Zone',
            'MHS': 'Residential - Mixed Housing Suburban Zone',
            'MHU': 'Residential - Mixed Housing Urban Zone',
            'THA': 'Residential - Terrace Housing and Apartment Building Zone'
               }
for postfix, zone in tqdm(postfix2name.items()):
    resid_gdf = resid[resid.ZONE_resol == zone].dissolve()
    parcels_sample[f'Hdist_{postfix}'] = parcels_sample.apply(lambda x: x.geometry.distance(resid_gdf.geometry[0]), axis=1)

CPU times: user 1min 42s, sys: 78.1 ms, total: 1min 42s
Wall time: 1min 42s


postfix = 'THA'
column = f'Hdist_{postfix}'
subsample = parcels_sample.sample(10)
subsample['buffer'] = subsample.apply(lambda x: x.geometry.buffer(x[column]), axis=1)
subsample['geometry'] = subsample['buffer']
ax = subsample.plot(color='red', alpha=0.4)
resid[resid.ZONE_resol == postfix2name[postfix]].plot(ax=ax)
ctx.add_basemap(ax, crs=2193)


LA = gpd.read_file('input/Modified_Community_Boards_SHP.zip').to_crs(2193)
LA.sample(3)


parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample = gpd.sjoin(parcels_sample, LA[['Local_Area', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'Local_Area': 'Local_Area_name'})


sa2 = gpd.read_file('NZ-SA/statistical-area-2-2020-generalised.gdb').to_crs(2193)
sa2.sample(3)


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample = gpd.sjoin(parcels_sample, sa2[['SA22020_V1_00_NAME', 'SA22020_V1_00', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'SA22020_V1_00_NAME': 'SA22018_name', 'SA22020_V1_00': 'SA22018_code'})

CPU times: user 8.29 s, sys: 3.88 ms, total: 8.29 s
Wall time: 8.29 s


au2013 = gpd.read_file('input/area-unit-2013.gdb.zip').to_crs(2193)
au2013.sample(3)


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample = gpd.sjoin(parcels_sample, au2013[['AU2013_V1_00_NAME', 'AU2013_V1_00', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'AU2013_V1_00_NAME': 'AU2013_name', 'AU2013_V1_00': 'AU2013_code'})

CPU times: user 6.96 s, sys: 0 ns, total: 6.96 s
Wall time: 6.96 s


mb2018 = gpd.read_file('input/meshblock-2018-clipped-generalised.gdb.zip').to_crs(4326)
mb2018.sample(3)


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid']
parcels_sample = gpd.sjoin(parcels_sample, mb2018[['MB2018_V1_00', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'MB2018_V1_00': 'MB2018_code'})

CPU times: user 11.3 s, sys: 8.2 ms, total: 11.3 s
Wall time: 11.3 s


mb2013 = gpd.read_file('input/meshblock-2013.gdb.zip').to_crs(4326)
mb2013.sample(3)


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample = gpd.sjoin(parcels_sample, mb2013[['MeshblockNumber', 'geometry']]).drop(columns=['index_right'])
parcels_sample = parcels_sample.rename(columns={'MeshblockNumber': 'MB2013_code'})

CPU times: user 11.4 s, sys: 7.05 ms, total: 11.4 s
Wall time: 11.4 s


parcels_sample['geometry'] = parcels_sample['geometry_centroid']
parcels_sample = parcels_sample.to_crs(2193)
parcels_sample.sample(3)


%%time
coastline = gpd.read_file('input/lds-nz-coastline-mean-high-water-FGDB.zip!nz-coastline-mean-high-water.gdb').to_crs(2193)
coastline = coastline.cx[bounds2193['x1']:bounds2193['x2'], bounds2193['y1']:bounds2193['y2']]
coastline_dissolved = coastline.dissolve()


%%time
parcels_sample['Hdist_coast'] = parcels_sample.apply(lambda x: x.geometry.distance(coastline_dissolved.geometry[0]), axis=1)

CPU times: user 14.3 s, sys: 0 ns, total: 14.3 s
Wall time: 14.3 s


# if distance work, then red circles should extend to the nearest coastline, and no further
subsample = parcels_sample.sample(10)
subsample['coast_buffer'] = subsample.apply(lambda x: x.geometry.buffer(x.Hdist_coast), axis=1)
subsample['geometry'] = subsample['coast_buffer']
ax = subsample.plot(color='red', alpha=0.4)
coastline.cx[1.7e6:1.8e6, 5.85e6:5.97e6].plot(ax=ax)

<AxesSubplot:>


roads = gpd.read_file('input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb').to_crs(2193)
roads = roads.cx[bounds2193['x1']:bounds2193['x2'], bounds2193['y1']:bounds2193['y2']]
highways = roads[~roads.hway_num.isna()]
highways_dissolved = highways.dissolve()


ax = highways.plot()
ctx.add_basemap(ax, crs=2193)


parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample['Hdist_motorway'] = parcels_sample.apply(lambda x: x.geometry.distance(highways_dissolved.geometry[0]), axis=1)


# if distance work, then red circles should extend to the nearest coastline, and no further
subsample = parcels_sample.sample(10)
subsample['highway_buffer'] = subsample.apply(lambda x: x.geometry.buffer(x.Hdist_motorway), axis=1)
subsample['geometry'] = subsample['highway_buffer']
ax = subsample.plot(color='red', alpha=0.4)
highways.plot(ax=ax)

<AxesSubplot:>


railroads = gpd.read_file('input/lds-nz-railway-centrelines-topo-150k-SHP.zip').to_crs(2193)
railroads = railroads.cx[bounds2193['x1']:bounds2193['x2'], bounds2193['y1']:bounds2193['y2']]
railroads_dissolved = railroads.dissolve()


parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample['Hdist_rail'] = parcels_sample.apply(lambda x: x.geometry.distance(railroads_dissolved.geometry[0]), axis=1)


# if distance work, then red circles should extend to the nearest coastline, and no further
subsample = parcels_sample.sample(10)
subsample['rail_buffer'] = subsample.apply(lambda x: x.geometry.buffer(x.Hdist_rail), axis=1)
subsample['geometry'] = subsample['rail_buffer']
ax = subsample.plot(color='red', alpha=0.4)
railroads_dissolved.plot(ax=ax)
ctx.add_basemap(ax, crs=2193)


skytower = [-36.84838748948485, 174.7621736911587]
skytower = gpd.points_from_xy(x=[skytower[1]], y=[skytower[0]])
skytower = gpd.GeoDataFrame([{"name": "Skytower", "value": 1}], geometry=skytower, crs="EPSG:4326").to_crs(epsg=2193)


parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample['Hdist_skytower'] = parcels_sample.apply(lambda x: x.geometry.distance(skytower.geometry[0]), axis=1)


# if distance works, then red circles should extend to the nearest sky tower, and no further
subsample = parcels_sample.sample(10)
subsample['skytower_buffer'] = subsample.apply(lambda x: x.geometry.buffer(x.Hdist_skytower), axis=1)
subsample['geometry'] = subsample['skytower_buffer']
ax = subsample.plot(color='red', alpha=0.2)
skytower.plot(ax=ax, color='black')
coastline.cx[1.7e6:1.8e6, 5.85e6:5.97e6].plot(ax=ax)

<AxesSubplot:>


spha = gpd.read_file('input/AC_Special_Housing_Area.zip').to_crs(2193)
spha_dissolved = spha.dissolve()
assert(len(spha_dissolved) == 1)


%%time
parcels_sample['geometry'] = parcels_sample['geometry_centroid_2193']
parcels_sample['SpHA_indicator'] = parcels_sample.apply(lambda x: spha_dissolved.geometry.contains(x.geometry), axis=1)


# if distance works, then red circles should extend to the nearest sky tower, and no further
subsample = parcels_sample.sample(500)
ax=subsample.plot(column='SpHA_indicator', alpha=0.6)
plt.ylim((5.89e6, 5.95e6))
plt.xlim((1.73e6, 1.78e6))
spha_dissolved.boundary.plot(ax=ax)
ctx.add_basemap(ax, crs=spha_dissolved.crs)

	OBJECTID	MXLOCATION	designvolt	status	descriptio	type	Symbol	Shape__Len	GlobalID	geometry
182	193	HEN-OTA-A	220	COMMISSIONED	Henderson - Otahuhu A	TRANSLINE	220 TRANSLINE	33182.373990	f7f50eff-0bec-490a-ab7f-db29aeb435d4	MULTILINESTRING ((1762946.160 5911750.110, 176...
180	191	HAM-MER-A	110	COMMISSIONED	Hamilton - Meremere A	TRANSLINE	110 TRANSLINE	58844.489544	e815b7ec-612a-493d-b9fa-322b338574e1	MULTILINESTRING ((1802908.070 5816213.260, 180...
223	2701	OTA-WKM-B	220	COMMISSIONED	Otahuhu - Whakamaru B	TRANSLINE	220 TRANSLINE	188931.028113	71f4fde0-6452-4acf-b612-059b1aab0ea4	MULTILINESTRING ((1844421.478 5743965.291, 184...

	OBJECTID	MXLOCATION	designvolt	status	descriptio	type	Symbol	Shape__Len	GlobalID	geometry
3	4	ALB-HPI-A	220	COMMISSIONED	Albany - Huapai A	TRANSLINE	220 TRANSLINE	14887.805068	8cb81a31-2c9c-497c-a8fc-646614b9290d	MULTILINESTRING ((1737714.820 5930696.980, 173...
4	5	ALB-SVL-A	220	COMMISSIONED	Albany - Silverdale A	TRANSLINE	220 TRANSLINE	14346.263293	ba383664-b20a-4f21-ae5e-0116474be615	MULTILINESTRING ((1750911.000 5932791.000, 175...
9	10	ARI-HAM-A	110	COMMISSIONED	Arapuni - Hamilton A	TRANSLINE	110 TRANSLINE	47281.457509	21c1d657-5f78-466f-b37d-bd48ff5ceef6	LINESTRING (1831824.000 5782923.000, 1831615.2...
10	11	ARI-HAM-B	110	COMMISSIONED	Arapuni - Hamilton B	TRANSLINE	110 TRANSLINE	47307.999476	571fd289-0092-4620-87ca-59437d29c741	MULTILINESTRING ((1831813.100 5782888.991, 183...
30	32	BOB-MER-A	110	COMMISSIONED	Bombay - Meremere A	TRANSLINE	110 TRANSLINE	15964.930587	b16a0624-3f79-45d5-b408-c51b36c3c594	MULTILINESTRING ((1783299.000 5868181.000, 178...
...	...	...	...	...	...	...	...	...	...	...
199	210	OTA-WKM-C	220	COMMISSIONED	Otahuhu - Whakamaru C	TRANSLINE	220 TRANSLINE	192292.688522	9e1d67d0-6186-46a1-a1c3-f79c2dff1440	MULTILINESTRING ((1844420.000 5743894.000, 184...
210	222	HAM-DEV-A	220	COMMISSIONED	Hamilton - Deviation A	TRANSLINE	220 TRANSLINE	6198.347734	101c45bc-9dd4-4815-b147-f4efbf6438e4	MULTILINESTRING ((1802971.630 5816033.640, 180...
215	227	KPU-WKO-A	110	COMMISSIONED	Kopu - Waikino A	TRANSLINE	110 TRANSLINE	33561.365858	1e4bd1ef-9b89-4ce7-86f3-0e3befb7bd3f	MULTILINESTRING ((1845525.290 5855448.280, 184...
222	2700	OTA-WKM-A	220	COMMISSIONED	Otahuhu - Whakamaru A	TRANSLINE	220 TRANSLINE	188787.769633	7b5e1909-cca6-4f81-866e-d32ea6bbcc3d	MULTILINESTRING ((1844422.210 5744003.450, 184...
223	2701	OTA-WKM-B	220	COMMISSIONED	Otahuhu - Whakamaru B	TRANSLINE	220 TRANSLINE	188931.028113	71f4fde0-6452-4acf-b612-059b1aab0ea4	MULTILINESTRING ((1844421.478 5743965.291, 184...

	id	appellation	affected_surveys	parcel_intent	topology_type	statutory_actions	land_district	titles	survey_area	calc_area	geometry	geometry_centroid_4326	geometry_polygon_4326	geometry_centroid_2193	geometry_polygon_2193	LINZ_parcel_vertices_lon	LINZ_parcel_vertices_lat	index_right	OBJECTID	TYPE	TYPE_resol	SUBTYPE	SCHEDULE	NAME	VALIDATION	VALIDATI00	last_edi00	created_da	VERSIONSTA	VERSIONS00	DocumentUR	GlobalID	SHAPE_Leng	SHAPE_Area
332148	6611744	Section 15 SO 317214	SO 317214	Road	Primary	[Create] Road, Limited Access Road and State H...	North Auckland	None	594.0	600.0	MULTIPOLYGON (((1749971.336 5940738.679, 17499...	POINT (174.67774 -36.66742)	MULTIPOLYGON (((174.67817 -36.66715, 174.67822...	POINT (1749932.255 5940710.038)	MULTIPOLYGON (((1749971.336 5940738.679, 17499...	[174.678168; 174.678224; 174.677304; 174.67725...	[-36.667152; -36.667189; -36.667681; -36.66764...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
138662	4971220	Lot 71 DP 4276	DP 4276	DCDB	Primary	None	North Auckland	NA184/150	612.0	613.0	MULTIPOLYGON (((1755077.168 5915667.766, 17550...	POINT (174.74041 -36.89243)	MULTIPOLYGON (((174.74038 -36.89224, 174.74055...	POINT (1755079.119 5915646.345)	MULTIPOLYGON (((1755077.168 5915667.766, 17550...	[174.74038; 174.740546; 174.740433; 174.740267...	[-36.892238; -36.892272; -36.892623; -36.89258...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
276820	6555480	Lot 2 DP 204937	DP 204937	Fee Simple Title	Primary	None	North Auckland	NA133B/442	475.0	475.0	MULTIPOLYGON (((1771091.292 5912435.935, 17710...	POINT (174.92060 -36.91862)	MULTIPOLYGON (((174.92074 -36.91859, 174.92076...	POINT (1771078.208 5912432.399)	MULTIPOLYGON (((1771091.292 5912435.935, 17710...	[174.920744; 174.920757; 174.920695; 174.92073...	[-36.918588; -36.91862; -36.918719; -36.918737...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
31004	4744873	Lot 106 DP 57498	DP 57498	DCDB	Primary	None	North Auckland	NA11B/1184	657.0	657.0	MULTIPOLYGON (((1751708.537 5925603.069, 17517...	POINT (174.70078 -36.80340)	MULTIPOLYGON (((174.70061 -36.80326, 174.70078...	POINT (1751723.367 5925587.744)	MULTIPOLYGON (((1751708.537 5925603.069, 17517...	[174.700606; 174.70078; 174.700945; 174.700846...	[-36.803262; -36.803205; -36.803533; -36.80356...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
494255	4696478	Lot 49 DP 71916	DP 71916	DCDB	Primary	None	North Auckland	NA38A/689	675.0	675.0	MULTIPOLYGON (((1755689.311 5929951.217, 17557...	POINT (174.74443 -36.76359)	MULTIPOLYGON (((174.74432 -36.76344, 174.74465...	POINT (1755698.853 5929934.142)	MULTIPOLYGON (((1755689.311 5929951.217, 17557...	[174.744322; 174.744652; 174.744543; 174.74421...	[-36.763436; -36.763582; -36.763741; -36.76359...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
24891	4732360	Lot 1 DP 189392	DP 189392	DCDB	Primary	None	North Auckland	NA118B/374	860.0	793.0	MULTIPOLYGON (((1733103.033 5930453.056, 17331...	POINT (174.49143 -36.76224)	MULTIPOLYGON (((174.49127 -36.76235, 174.49127...	POINT (1733117.978 5930465.529)	MULTIPOLYGON (((1733103.033 5930453.056, 17331...	[174.491265; 174.491267; 174.491269; 174.49126...	[-36.762354; -36.762346; -36.762331; -36.76231...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
25636	4733918	Lot 105 DP 59567	DP 59567	DCDB	Primary	None	North Auckland	NA14C/462	703.0	703.0	MULTIPOLYGON (((1770070.630 5918154.599, 17700...	POINT (174.90802 -36.86744)	MULTIPOLYGON (((174.90801 -36.86725, 174.90821...	POINT (1770070.878 5918133.805)	MULTIPOLYGON (((1770070.630 5918154.599, 17700...	[174.908008; 174.908208; 174.908026; 174.90782...	[-36.867249; -36.867551; -36.867628; -36.86732...	5.0	6.0	3	Height Sensitive Areas	None	None	None	3	Valid and Public	20161109010	20160530235	4	Operative	None	{7EE3FD70-DD2B-4138-B975-6CA1301626A8}	9595.090952	1.303548e+06
25636	4733918	Lot 105 DP 59567	DP 59567	DCDB	Primary	None	North Auckland	NA14C/462	703.0	703.0	MULTIPOLYGON (((1770070.630 5918154.599, 17700...	POINT (174.90802 -36.86744)	MULTIPOLYGON (((174.90801 -36.86725, 174.90821...	POINT (1770070.878 5918133.805)	MULTIPOLYGON (((1770070.630 5918154.599, 17700...	[174.908008; 174.908208; 174.908026; 174.90782...	[-36.867249; -36.867551; -36.867628; -36.86732...	6.0	7.0	4	Viewshafts	None	B6	Browns Island	3	Valid and Public	20161109010	20160531002	4	Operative	None	{4921E0E2-AFE9-4BAA-9D17-CC5669FC65C1}	31274.305310	3.767934e+07
369084	7503911	Part Lot 33 DP 19195	None	Fee Simple Title	Primary	None	North Auckland	NA46C/131	556.0	7.0	MULTIPOLYGON (((1753477.305 5920162.246, 17534...	POINT (174.72155 -36.85201)	MULTIPOLYGON (((174.72152 -36.85200, 174.72156...	POINT (1753479.593 5920160.987)	MULTIPOLYGON (((1753477.305 5920162.246, 17534...	[174.721524; 174.721557; 174.721568; 174.721524]	[-36.852002; -36.851997; -36.85204; -36.852002]	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
334475	4769767	Part Allot N35 PSH OF Ahuroa	SO 27	DCDB	Primary	None	North Auckland	NA14D/470	730609.0	120330.0	MULTIPOLYGON (((1743070.851 5963605.949, 17430...	POINT (174.59335 -36.46327)	MULTIPOLYGON (((174.59673 -36.46215, 174.59668...	POINT (1742765.292 5963486.442)	MULTIPOLYGON (((1743070.851 5963605.949, 17430...	[174.596732; 174.596678; 174.594612; 174.58897...	[-36.462145; -36.464002; -36.464061; -36.46422...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	WORKFLOWJO	last_edi00	CONTOUR	created_da	DocumentUR	GlobalID	GROUPZONE	GROUPZONE_	ID	NAME	OBJECTID	PARCEL_BAS	PARCEL_B00	PRECINCT	PRECINCT_r	PRECINCTGR	PRECINCT00	SCHEDULE	SHAPE_Area	SHAPE_Leng	STATUS	SUBPRECINC	SUBPRECI00	SUBTYPE	SUBTYPE_re	TYPE	TYPE_resol	VALIDATION	VALIDATI00	VERSIONSTA	VERSIONS00	ZONE	ZONE_resol	ZONEHEIGHT	geometry
58724	NaN	20161111011	None	20160718211	None	{0CCC66B0-9A15-477A-BE4F-FAE9C5FA9C15}	5	Coastal	None	None	58725	None	None	None	None	None	None	None	62.153950	58.346375	None	None	None	None	None	None	None	3	Valid and Public	4	Operative	59	Coastal - Coastal Transition Zone	NaN	POLYGON Z ((1756089.130 5973206.276 0.000, 175...
103279	NaN	20161111012	None	20160718211	None	{111B7716-6401-457F-A260-9BA18D7E0DB9}	5	Coastal	None	None	103280	None	None	None	None	None	None	None	372.662968	106.669373	None	None	None	None	None	None	None	3	Valid and Public	4	Operative	30	Coastal - General Coastal Marine Zone	NaN	POLYGON Z ((1810661.036 5993908.378 0.000, 181...
54081	NaN	20161111011	None	20160718211	None	{CD16D435-B3A5-4C26-BFC7-E6AF32E60231}	7	General	None	None	54082	None	None	None	None	None	None	None	41502.960080	1179.183755	None	None	None	None	None	None	None	3	Valid and Public	4	Operative	26	Strategic Transport Corridor Zone	NaN	POLYGON Z ((1737071.120 5973948.034 0.000, 173...

	WORKFLOWJO	last_edi00	CONTOUR	created_da	DocumentUR	GlobalID	GROUPZONE	GROUPZONE_	ID	NAME	OBJECTID	PARCEL_BAS	PARCEL_B00	PRECINCT	PRECINCT_r	PRECINCTGR	PRECINCT00	SCHEDULE	SHAPE_Area	SHAPE_Leng	STATUS	SUBPRECINC	SUBPRECI00	SUBTYPE	SUBTYPE_re	TYPE	TYPE_resol	VALIDATION	VALIDATI00	VERSIONSTA	VERSIONS00	ZONE	ZONE_resol	ZONEHEIGHT	geometry
20853	NaN	20161111011	None	20160718211	None	{2D83F680-9587-4D57-AF0B-CB7CA27C3D2C}	6	Special purpose zone	None	None	20854	None	None	None	None	None	None	None	22.286945	72.629617	None	None	None	None	None	None	None	3	Valid and Public	4	Operative	58	None	NaN	POLYGON Z ((1767684.860 5903714.023 33.721, 17...
121765	NaN	20161111010	None	20160718211	None	{FD3E8BB4-1979-42B5-9A77-5D04AE5190AF}	6	Special purpose zone	None	None	121766	None	None	None	None	None	None	None	4.114962	86.135341	None	None	None	None	None	None	None	3	Valid and Public	4	Operative	58	None	NaN	POLYGON Z ((1767684.083 5903713.690 33.721, 17...

	OBJECTID	Local_Area	geometry
8	10.0	Hibiscus Coast	MULTIPOLYGON (((1752023.352 5954803.281, 17520...
17	21.0	Papatoetoe	POLYGON ((1765105.421 5908611.893, 1765113.661...
2	25.0	Rodney-Kumeu-Riverhead	POLYGON ((1742532.808 5931237.574, 1742490.377...

	AU2013_V1_00	AU2013_V1_00_NAME	AREA_SQ_KM	LAND_AREA_SQ_KM	Shape_Length	geometry
1481	573903	Newlands North	0.804087	0.804087	5958.640337	MULTIPOLYGON (((174.82896 -41.21949, 174.83011...
958	525202	Kawakawa-Orere	105.216001	105.216001	57437.556145	MULTIPOLYGON (((175.14259 -36.93214, 175.14265...
608	597102	Inland Water-Lake Ellesmere South	63.304671	0.000000	75233.330853	MULTIPOLYGON (((172.57398 -43.76779, 172.57401...

	MB2018_V1_00	LANDWATER	LANDWATER_NAME	LAND_AREA_SQ_KM	AREA_SQ_KM	SHAPE_Length	geometry
52865	4011926	12	Mainland	0.173710	0.173710	2129.422679	MULTIPOLYGON (((174.91945 -36.94519, 174.92053...
23832	1429800	12	Mainland	0.021918	0.021918	731.394363	MULTIPOLYGON (((176.91682 -39.49291, 176.91702...
15554	0759520	12	Mainland	0.020695	0.020695	928.355735	MULTIPOLYGON (((174.91840 -37.01411, 174.91871...

	id	appellation	affected_surveys	parcel_intent	topology_type	statutory_actions	land_district	titles	survey_area	calc_area	geometry	geometry_centroid	geometry_polygon	SA22018_name	SA22018_code	MB2013_code	MB2018_code	AU2013_name	AU2013_code
102782	4896031	Lot 1 DP 61837	DP 61837	DCDB	Primary	None	North Auckland	NA17C/935	685.0	685.0	POINT (1773608.374 5914596.096)	POINT (174.94849 -36.89867)	MULTIPOLYGON (((174.94827 -36.89874, 174.94841...	Cockle Bay	153400	0655102	0655102	Cockle Bay	521502
119786	4931702	Lot 4 DP 166054	DP 166054	DCDB	Primary	None	North Auckland	NA100D/34	52130.0	50575.0	POINT (1783027.533 5889659.615)	POINT (175.06020 -37.12153)	MULTIPOLYGON (((175.06295 -37.12105, 175.06293...	Ararimu	166400	0815302	0815302	Hunua	521132
362342	7415715	Lot 23 DP 455616	DP 455616, DP 462513	Fee Simple Title	Primary	None	North Auckland	586776	103.0	103.0	POINT (1763989.858 5916157.717)	POINT (174.84025 -36.88632)	MULTIPOLYGON (((174.84016 -36.88629, 174.84020...	Stonefields West	144900	0465305	4000293	Stonefields	517202

Setup¶

1. LINZ parcel information¶

a. LINZ_parcel_ID Parcel identifier - “certificate of title”.¶

b. LINZ_parcel_centroid_lon¶

c. LINZ_parcel_centroid_lat¶

d. LINZ_parcel_vertices_lon vector of longitudes of the vertices of the matched LINZ parcels¶

e. LINZ_parcel_vertices_lat¶

f. LINZ_parcel_roadvertices_lon subvector of longitudes of parcel that sits adjacent to a road¶

g. LINZ_parcel_roadvertices_lat¶

h. LINZ_parcel_sides_zones AUP Zone Code of adjoining parcels (this includes residential, business, and rural zones; it should also include roads, water and open spaces)¶

i. i. LINZ_TRNSPWR_ohead_indicator Indicator (0 or 1) for LINZ parcel located under overhead transpower line. Leave blank otherwise. Note that ‘TRANSLINE’ denotes overhead transmission lines in the GIS dataset, while ‘CABLE’ denotes underground transmission lines and can be ignored.¶

i. ii. LINZ_TRNSPWR_ohead_name¶

j. i. LINZ_VWSHFT_ohead_indicator Indicator (0 or 1) for LINZ parcel located under viewshafts. Leave blank otherwise.¶

j. ii. LINZ_VWSHFT_ohead_name Name of the volcanic cone (e.g. Mt Albert). Leave blank if no viewshaft applies.¶

j. iii. LINZ_VWSHFT_ohead_ID OBJECTID of the viewshaft. Leave blank if no viewshaft applies.¶

2. AUP shapefile information. Calculations are based on centroid of matched parcels:¶

a. LINZmatch_AUP_name AUP zone (string) consent is located in (from the AUP shapefiles)¶

b. LINZmatch_AUP_code AUP Zone Code (numerical) consent is located in¶

c. Hdist_rural Minimum Haversine distance to nearest Rural Zone, plus¶

c. i. Hdist_rural_name Name of rural zone identified¶

c. ii. Hdist_rural_code Numeric code of rural zone identified¶

d. Hdist_bus Minimum Haversine distance to nearest Business Zone¶

d. i. Hdist_bus_name¶

d. ii. Hdist_bus_code¶

e. Hdist_resid Minimum Haversine distance to nearest Residential Zone,¶

e. i. Hdist_resid_name¶

e. ii. Hdist_resid_code¶

f. Hdist_SH Minimum Haversine distance to Residential - Single House Zone¶

g. Hdist_MHS Minimum Haversine distance to Residential - Mixed Housing Suburban Zone¶

h. Hdist_MHU Minimum Haversine distance to Residential - Mixed Housing Urban Zone¶

i. Zone Hdist_THA Minimum Haversine distance to Residential - Terrace Housing and Apartments¶

3. Local Area (LA) Information (information on the LA that the consent is located in):¶

a. Name (just the string is fine; I did not generate numeric codes for the local areas generated) Local_Area_name¶

4. 2018 Statistical Area 2 Information (information on the 2018SA2 that the consent is located in)¶

a. Name (string) SA22018_name¶

b. Numeric code SA22018_code¶

5. Area Unit Information (information on the 2013AU that the consent is located in):¶

a. Name (string) AU2013_name¶

b. Numeric code AU2013_code¶

6. 2018 Meshblock Information (information on the 2018MB that the consent is located in):¶

a. Code MB2018_code¶

7. 2013 Meshblock Information (information on the 2013MB that the consent is located in):¶

a. Code MB2013_code¶

8. Additional distance information from consent location¶

a. Minimum Haversine distance to coastline Hdist_coast¶

b. Minimum Haversinedistance to Open Space(set of AUP Zones) Hdist_open¶

c. Minimum Haversine distance to motorway Hdist_motorway¶

d. Minimum Haversine distance to main arterial road Hdist_main_road¶

e. Minimum Haversine distance to rail line Hdist_rail¶

f. Haversine distance to downtown (use Skytower coordinates) Hdist_skytower¶

9. Special Housing Area (SpHAs)¶

10. Ranged Address Indicator¶

11. Matching Type Indicators¶

a. LINZ_MATCH_CODE¶

b. LINZ_2ND_MATCH_CODE¶